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Abstract We use the planar circular restricted three-body 
problem in order to numerically investigate the orbital dy¬ 
namics of orbits of a spacecraft, or a comet, or an asteroid 
in the Saturn-Titan system in a scattering region around the 
Titan. The orbits can escape through the necks around the 
Lagrangian points L\ and L 2 or collide with the surface of 
the Titan. We explore all the four possible Hill’s regions de¬ 
pending on the value of the Jacobi constant. We conduct a 
thorough numerical analysis on the phase space mixing by 
classifying initial conditions of orbits and distinguishing be¬ 
tween three types of motion: (i) bounded, (ii) escaping and 
(iii) collisional. In particular, we locate the different basins 
and we relate them with the corresponding spatial distribu¬ 
tions of the escape and crash times. Our results reveal the 
high complexity of this planetary system. Furthermore, the 
numerical analysis shows a strong dependence of the prop¬ 
erties of the considered basins with the total orbital energy, 
with a remarkable presence of fractal basin boundaries along 
all the regimes. We hope our contribution to be useful in both 
space mission design and study of planetary systems. 

Keywords Restricted three body-problem; Escape dynam¬ 
ics 


1 Introduction 

Over the years many research works have been devoted on 
the issue of escaping or leaking particles from several types 
dynamical systems. In particular, the problem of escape is 
a classical problem in Hamiltonian nonlinear systems (e.g., 
Aguirre et al., 2001; Aguirre & Sanjuan, 2003; Aguirre et 
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al., 2009; Barrio et al., 2008; Blesa et al., 2012; Zotos, 2014) 
as well as in dynamical astronomy (e.g.. Hut & Bahcall, 
1983;Bleher etal., 1988;Benetet al., 1996,1998; de Moura 
& Letelier, 2000; Zotos, 2012), while escaping orbits in the 
classical Restricted Three-Body Problem (RTBP) is another 
typical example (e.g., Nagler, 2004, 2005; Davis & Howell, 
2011, 2012; de Assis & Terra, 2014; Zotos, 2015b). 

Nevertheless, the issue of escaping orbits in Hamiltonian 
systems is by far less explored than the closely related prob¬ 
lem of chaotic scattering (Bleher et al., 1988). In this sit¬ 
uation, a test particle coming from infinity approaches and 
then scatters off a complex potential. In a recently published 
paper (Zotos, 2015b) we investigated the orbital dynamics 
in the Copenhagen problem where one of the two primary 
bodied was a oblate spheroid. There the scattering region on 
the configuration (x, y) space was a rectangular region ex¬ 
tending around both primaries and we managed to classified 
initial conditions of orbits into different categories (i.e., or¬ 
bits that crash into one of the primaries, escaping orbits and 
bounded orbits of several types). In the current work on the 
other hand, the scattering region of initial conditions of or¬ 
bits is confined only around the Titan and especially between 
the Lagrangian points Li and L 2 . 

Escaping and trapped motion of stars in stellar systems 
is an another issue of great importance. In a previous article 
(Zotos, 2012), we explored the nature of orbits of stars in a 
galactic-type potential, which can be considered to describe 
local motion in the meridional (R, z) plane near the central 
parts of an axially symmetric galaxy. It was observed, that 
apart from the trapped orbits there are two types of escap¬ 
ing orbits, those which escape fast and those which need 
to spend vast time intervals inside the equipotential surface 
before they find the exit and eventually escape. The escape 
dynamics and the dissolution process of a star cluster em¬ 
bedded in the tidal field of a parent galaxy was investigated 
in Ernst et al. (2008). Conducting a thorough scanning of 
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the available phase space the authors managed to obtain the 
basins of escape and the respective escape rates of the orbits, 
revealing that the higher escape times correspond to initial 
conditions of orbits near the fractal basin boundaries. The 
investigation was expanded into three dimensions in Zotos 
(2015a) where we revealed the escape mechanism of three- 
dimensional orbits in a tidally limited star cluster. Further¬ 
more, Ernst & Peters (2014) explored the escape dynamics 
in the close vicinity of and within the critical area in a two- 
dimensional barred galaxy potential, identifying the escape 
basins both in the phase and the configuration space. The 
numerical approach of the above-mentioned papers served 
as the basis of this work. 

Many times in the past the RTBP has played an essential 
role in several areas of celestial mechanics and dynamical 
astronomy. For instance, the modern applications to space 
flight missions (Koon et al., 2000b; Gomez et al., 2001a,b,c,d, 
2004; Stuchi et al., 2008) are more numerous than the the¬ 
oretical classical applications. Moreover, today many issues 
in space dynamics are of paramount importance and interest. 
The applications of RTBP include the launching of artificial 
satellites in the solar system and they also form the basis of 
several planetary theories. 

Titan, Saturn’s largest moon, is undoubtedly a treasure 
trove of potential space discovery. Having an environment 
which is comparable to that of the primordial Earth and an 
atmosphere denser than that of any other moon in the so¬ 
lar system makes it a perfect candidate for space missions. 
In particular,the Cassini-Huygens mission’s recent discov¬ 
ery of hydrocarbon lakes on the Titan’s surface (Koon et al., 
2000a) has so inspired the National Aeronautics and Space 
Administration (NASA) that they plan to deploy a space 
mission to Titan once more in roughly one decade from now. 
In the present contribution we shall explore the orbital dy¬ 
namics of the Saturn-Titan system applying as a guide the 
numerical approach used in de Assis & Terra (2014) (here¬ 
after Paper I) for the Earth-Moon system. 

The structure of the paper is as follows: In Section 2 we 
present in detail the properties of the dynamical model. All 
the fundamental concepts and the computational methods 
we used in order to determine the classification of the orbits 
are described in Section 3. In the following Section, we con¬ 
duct a thorough and systematic numerical investigation re¬ 
vealing the overall orbital structure of the Saturn-Titan plan¬ 
etary system. Our paper ends with Section 5 where the dis¬ 
cussion and the main conclusions of this work are given. 


2 The dynamical model 

Let us briefiy recall the basic properties of the planar cir¬ 
cular restricted threebody problem (PCRTBP) (Szebehely, 
1967). The two main bodies, called primaries Pi and P 2 


move on circular orbits around their common center of grav¬ 
ity. The third body (test particle) moves in the same plane 
under the gravitational field of the two primaries. It is as¬ 
sumed that the motion of the two primaries is not perturbed 
by the third body since the third body’s mass is much smaller 
with respect to the masses of the two primaries. The non- 
dimensional masses of the two primaries are 1 - yu and yu, 
where yu = m 2 l(mi m 2 ) is the mass ratio and mi > m 2 . For 
the Saturn-Titan (ST) system, fi = 0.0002461294 (Sharma 
& Subba Rao, 1976; Beevi & Sharma, 2012). We choose as 
a synodic barycenter reference frame a rotating coordinate 
system where the origin is at the center of mass of Pi - P 2 , 
while the centers Ci and C 2 of the two primaries are located 
at (-yu, 0) and (1 - yu, 0), respectively. 

The total time-independent effective potential function 
in the rotating frame is 


Q(x,y) = 


ri r2 2rj 2 ' ' 


where 


( 1 ) 


ri = 





ri 



+ ix-lf+y^. 


-r 


3Ai 


( 2 ) 


are the distances to the respective primaries and the angular 
velocity (n). 

The larger primary body Pi (Saturn) is an oblate spheroid 
(Beatty et al., 1999) and Ai is the oblateness coefficient which 
is defined as 

_ {REf-{RPf 

‘ 5^2 ’ ( 3 ) 

where RE and RP are the equatorial and polar radius, respec¬ 
tively of the oblate primary, while R is the distance between 
the centers of the two primaries. According to Table 1 of 
Sharma & Subba Rao (1976) the oblateness coefficient for 
the Saturn-Titan system is Ai = 0.0000981153. The study 
of oblateness includes the series of works of Markellos et 
al. (2002); Kalantonis et al. (2005,2006, 2008); Kalvouridis 
& Gousidou-Koutita (2012); Perdiou et al. (2012); Sharma 
& Subba Rao (1979, 1986); Subba Rao & Sharma (1988, 
1997); Sharma (1981,1987, 1989, 1990) by considering the 
more massive primary as an oblate spheroid with its equa¬ 
torial plane co-incident with the plane of motion of the pri¬ 
maries. 

The scaled equations of motion describing the motion of 
the third body in the corotating frame read (Beevi & Sharma, 
2012) 


dQ ^ dQ 
X - 2ny = —, y 2nx = —. 

ox oy 


( 4 ) 
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Fig. 1 The isolines contours of the constant potential and the location 
of the Lagrangian points around smaller primary (Titan). The interior 
region around Titan is indicated in white, the region around Saturn is 
colored in green, the exterior region is shown in red, while the forbid¬ 
den regions of motion are marked with grey. 


The dynamical system (4) admits the well know Jacobi in¬ 
tegral 

J{x, y, X, y) = 2Q{x, y) - (^x^ -\-y^^ = C, (5) 

where x and y are the velocities, while C is the Jacobi con¬ 
stant which is conserved and defines a three-dimensional in¬ 
variant manifold in the total four-dimensional phase space. 
Therefore, an orbit with a given value of it’s energy integral 
is restricted in its motion to regions in which C < 2Q(x,y), 
while all other regions are forbidden to the third body. The 
energy value E is related with the Jacobi constant by C = 
-2E. 

The PCRTBP with oblateness has five equilibria known 
as Lagrangian points at which 


dQ _ dQ 
dx dy 


( 6 ) 


Three of them, Li, L2, and L3, are collinear points located in 
the v-axis, while the other two L 4 and L 5 are called triangular 
points and they are located in the vertices of an equilateral 
triangle. The central stationary point Li is a local minimum 
of the potential function Q{x,y), while the stationary points 
L 2 and L 3 are saddle point. Li is between the primaries, L 2 is 
at the right side of P 2 , while L 3 is at the left side of Pi. The 
points L 4 and L 5 on the other hand, are local maxima of the 
potential function, enclosed by the banana-shaped isolines. 
The projection of the four-dimensional phase space onto the 
configuration (or position) space (v,y) is called the Hill’s 
regions. The boundaries of these Hill’s regions are called 


Zero Velocity Curves (ZVCs) because they are the locus in 
the configuration (v, y) space where the kinetic energy van¬ 
ishes. The Jacobi constant values at are denoted by C^. 
For the Saturn-Titan system we have: Ci = 3.0164219315, 
C2 = 3.0160936366, C3 = 3.0004914646, and C4 = C5 = 
2.9999991590. In Fig. 1 we see that the presence of forbid¬ 
den regions in the PCRTBP permits the definition of four 
subsets of the configuration (v, y) plane when C 2 < C < Cy. 
the interior region around Titan (white), the region around 
Saturn (green), the exterior region (red), and the forbidden 
regions (gray). 


3 Computational methods 

For exploring the orbital dynamics in the Saturn-Titan sys¬ 
tem, we need to define samples of initial conditions of orbits 
whose properties will be identified. For this purpose we con¬ 
sider dense, uniform grids of 1024 x 1024 initial conditions 
(•^0, To) regularly distributed on the configuration (x, y) plane 
inside the area allowed by the Jacobi constant C. Following 
a typical approach, all orbits are launched with initial condi¬ 
tions inside a certain region, called scattering region, which 
in our case is xl^ < x < xl 2 and -0.1 < y < 0.1. For all 
orbits Vo = 0 , while the value of y'o is always obtained from 
the Jacobi integral (5) as y'o = y(vo, yo, C) > 0. 

In the PCRTBP system the configuration space extends 
to infinity thus making the identification of the type of mo¬ 
tion of the third body for specific initial conditions a rather 
demanding task. Depending on the Jacobi constant the re¬ 
gion around the Titan can be connected or not to the Sat¬ 
urn realm through the neck around Li or to the exterior 
realm through the neck around L 2 . Therefore initial con¬ 
ditions of orbits around the Titan can be categorized into 
basins \ namely: (i) the bounded basin containing orbits which 
remain in the Titan realm, (ii) the Saturn realm basin consti¬ 
tuted by orbits which escape through Li, (iii) the exterior 
realm basin corresponding to orbits that escape through L 2 , 
(iv) the collisional basin or orbits which collide with the Ti¬ 
tan. 

Now we need to define appropriate numerical criteria for 
distinguishing between these four types of motion. In order 
to consider a more realistic approach, we assume that the 
Titan is a finite body, taking into account its mean radius 
approximately by 2576 km (about 6.7 x 10“^ dimensionless 
length units). Therefore, if an orbit reaches the surface of the 
Titan, its numerical integration ends thus producing an orbit 
leaking in the configuration space. Furthermore, an escaping 
orbit to the Saturn realm must satisfy the conditions v < 
vli - di, with Si = 0.07 and the third body inside he circle 
around the Saturn of radius ri < Ix^ - v^J. In the same 

^ The set of initial conditions of orbits which lead to a certain final 
state (escape, collision or bounded motion) is defined as a basin. 
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vein, an escaping orbit to the exterior realm must fulfill the 
conditions x > + S 2 , with S 2 = 0.05, or |}7| > or, if 

X < xl^-Si, ri > |xl 3 -xpj I (i.e., the third body is outside the 
circle just defined). Here we must clarify that the tolerances 
and S 2 were included in the escape criteria in an attempt 
to avoid that the unstable Lyapunov orbits are incorrectly 
classified as escaping orbits. 

In our numerical integrations the maximum time tf em¬ 
ployed is 5000 dtu (dimensionless time units), correspond¬ 
ing to about 64.315 yr. Usually, the vast majority of orbits 
need considerable less time to escape from the system (ob¬ 
viously, the numerical integration is effectively ended when 
an orbit moves outside the system’s disk and escapes). Nev¬ 
ertheless, we decided to use such a long integration time just 
to be sure that all orbits have enough time in order to escape. 
Remember, that there are the so called “sticky orbits” which 
behave as regular ones during long periods of time. Here we 
should clarify, that orbits which do not escape or collide to 
Titan after a numerical integration of 5000 dtu are consid¬ 
ered as bounded regular orbits. 

The equations of motion (4) for the initial conditions 
of all orbits are forwarded integrated using a double pre¬ 
cision Bulirsch-Stoer FORTRAN 77 algorithm (e.g.. Press et 
al. (1992)) with a small time step of order of 10“^, which 
is sufficient enough for the desired accuracy of our compu¬ 
tations. Here we should emphasize, that our previous nu¬ 
merical experience suggests that the Bulirsch-Stoer integra¬ 
tor is both faster and more accurate than a double precision 
Runge-Kutta-Fehlberg algorithm of order 7 with Cash-Karp 
coefficients. Throughout all our computations, the Jacobi in¬ 
tegral (Eq. (5)) was conserved better than one part in 10“^^ 
although for most orbits it was better than one part in 10“^^. 
For collisional orbits where the test body moves inside a re¬ 
gion of radius 10“^ around the Titan the Lemaitre’s global 
regularization method is applied. 

4 Numerical results 

The value of the Jacobi constant defines four different cases 
regarding the structure of the Hill’s regions. In particular 

- Case I: C > C\ All the transport channels around Li 
and L 2 are closed, so there is only the bounded and the 
collisional basins. 

- Case II: Cl > C > C 2 Only the neck around Li is open. 
Therefore, the Saturn realm basin joins the two basins of 
the previous case. 

- Case III: C 2 > C > C 3 Both necks around Li and L 2 are 
open. Now we have the exterior realm basin along with 
the three basins of the second case. 

- Case IV: C < C3 All four basins are present. In this case 
however, the forbidden regions are reduced (C3 > C > 
C4) or completely absent (C < C4). 


4.1 Case I: C > Ci 

In Fig. 2 we present the basin diagrams for four values of 
the Jacobi constant. The outermost black solid line is the 
ZVC which is defined as Q(x,y) = Cjl. Examining these 
diagrams we see that in Case I there are only bounded and 
collisional orbits. In particular, for C = 3.04, that is an en¬ 
ergy level very close to the Titan’s surface, all the available 
configuration space is covered entirely by initial conditions 
which correspond to collisional orbits^. When C = 3.03 two 
stability islands emerge inside the vast collisional basin. The 
smaller stability island, located on the left side of the Titan, 
surrounds a periodic orbit around P 2 which is symmetric to 
a reflection over the x-axis (y = 0 ) and is traveled in counter¬ 
clockwise sense, hence, prograde with respect to the rotating 
coordinate system. The larger stability island on the other 
hand, situated on the right side of the Titan, surrounds a 
retrograde (clockwise) symmetric periodic orbit around P 2 . 
Simo & Stuchi (2000) proved for the planar Hill problem 
(see also de Assis & Terra (2014)) that the stability island 
on the left side of the Titan is much more stable than that on 
the right side in relation to C variation. As the value of the 
Jacobi constant decreases we see in Fig. 2c that for C = 3.02 
both stability islands grow in size. For the lowest value of C, 
for instance C = 3.01643 (Fig. 2d), we observe that the size 
of the stability island located on the left side of the Titan re¬ 
mains almost unperturbed, while the stability region on the 
right side considerably reduces. 


4.2 Case II: Ci > C > C2 

When C < Ci the neck around Li opens thus allowing an 
orbit to enter the Saturn realm. Indeed for C = 3.01641 we 
see in Fig. 3a that there is an escape channel in the left side 
of the ZVC. However, even though that the escape is possi¬ 
ble we did not detected any initial conditions of orbits that 
escape through Li and therefore, the configuration plane is 
still composed of only bounded and collisional initial condi¬ 
tions. Perhaps for such values of C the width of the escape 
channel is too small, so all orbits collide to the Titan’s sur¬ 
face before they manage to escape. For C = 3.01633 the 
width of the escape channel doubles in size and one can ob¬ 
serve the presence of thin filaments of escape basins around 
the stability island located in the right side of the Titan. As 
the value of C decreases the width of the escape channel in¬ 
creases thus allowing more and more orbits to escape to the 
Saturn realm. Apart from the increase to the escape basins 
we observe that the stability island in the right side of the 
Titan reduces in size. At the lowest value of the Jacobi con¬ 
stant in Case II, that is C = 3.01610, we identify in Fig. 3d 

^ As in Paper I we also defined initial conditions inside the radius of 
the Titan which are obviously collisional orbits. 
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Fig. 2 Basin diagrams for Case L (a-upper left): C = 3.04, (b-upper right): C = 3.03, (c-lower left): C = 3.02 and (d-lower right): C = 3.01643. 
The color code is as follows: bounded basin (gray), collisional basin (cyan). 





Fig. 3 Basin diagrams for Case II. (a-upper left): C = 3.01641, (b-upper right): C = 3.01633, (c-lower left): C = 3.01625 and (d-lower right): 
C = 3.01610. The color code is as follows: bounded basin (gray), collisional basin (cyan), Saturn real basin (green). 

Here it should be emphasized that the numerical inte¬ 
gration is stopped as soon as the third body passes Li thus 
entering the Saturn realm. However, if we do not stop the 
numerical integration some orbits that initially entered the 


well defined escape basins inside the vast collisional basin, 
while on the right side of the diagram, only a minor bounded 
basin is present. 
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Fig. 4 Basin diagrams for Case III. (a): C = 3.01607, (b): C = 3.0159, (c): C = 3.0154, (d): C = 3.013, (e): C = 3.007, and (f): C = 3.001. The 
color code is as follows: bounded basin (gray), collisional basin (cyan), Saturn realm basin (green), exterior realm basin (red). 


Saturn realm may return to the Titan region, or crash into 
on of the primary bodies, or even enter the exterior region 
and then escape to infinity. In our calculations we follow the 
approach used in de Assis & Terra (2014) and we consider 
an orbit to escape to Saturn realm if the test body passes Li 
even if its true asymptotic behaviour at very long time limit 
is different. 

4.3 Case III: C 2 > C > C 3 

This case constitutes the Hill’s regions with the most inter¬ 
est in point of view of planetary systems. Furthermore, it 
has many practical applications, such as orbit determination 
of spacecraft mainly based on many-models and also the 
phenomenon of temporary capture of comets or asteroids 
around planets in our solar system. The exterior realm rises 
when C < C 2 . In Fig. 4a we see that for C = 3.01607 both 


necks around Li and L 2 are open. Once more, the vast ma¬ 
jority of the configuration (v,y) space is covered by initial 
conditions of orbits which collide to the Titan. It should be 
noticed that a small escape basin corresponding to the exte¬ 
rior realm is located at the right outer part of the plane very 
close to L 2 . When C = 3.0159 it is seen in Fig. 4b that ad¬ 
ditional exterior realm basins emerge inside the collisional 
basin, while the stability island in the right side of the Ti¬ 
tan is no longer present. As we proceed to smaller values of 
the Jacobi constat (see Figs. 4(c-e)) the extent of the exterior 
realm basins grows significantly, while the collisional basin 
shrinks thus confined mainly near the central region of the 
plane. In particular we see in Fig. 4e that when C = 3.007 
about 45% of the available configuration plane is occupied 
by initial conditions of orbits that escape through L 2 to the 
exterior region. Finally when C = 3.001 it is evident that 
the exterior realm basin take over the configuration plane. 
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Fig. 5 Basin diagrams for Case IV. (a-upper left): C = 3, (b-upper right): C = 2.999, (c-lower left): C = 2.995 and (d-lower right): C = 2.98. The 
color code is the same as in Fig. 4. 

leading to the predominance of this type of orbits at the end plane is occupied by initial conditions composing the uni- 
of Case III. Moreover, the extent of the collisional and the fled exterior realm basin. We note that the collisional basin 
Saturn realm basins is considerably reduced, while a small has a stream-like shape which is initiated at the center of the 
stability island in the left side of the Titan survives. plane around the radius of the Titan. 


4.4 Case IV: C < C 3 

In this last case, there are two subintervals delimited by two 
critical values of the jacobi constant. In particular, at C = C 3 
the neck associated with the Lagrangian point L 3 opens. 
Therefore, in the first subcase, C 3 > C > C 4 the ZVC and 
the forbidden regions are still present on the configuration 
space, while for C < C 4 the forbidden regions disappears 
and all the (x, y) plane is available for motion. Fig. 5a shows 
the structure of the configuration plane for C = 3. We ob¬ 
serve that about 80% of the plane is covered by initial con¬ 
ditions of orbits corresponding to escape through L 2 . With 
a much closer look we can identify a small stability island 
of bounded motion in the left side of the Titan. It should 
be pointed out that the vicinity of this stability region is 
highly fractal. This stability island is no longer present in 
Fig. 5b where C = 2.999 thus the configuration plane con¬ 
sists entirely of escape and collisional basins. As we pro¬ 
ceed to lower values of the Jacobi constant the extent of the 
Saturn realm basin is reduced, while that of the collisional 
basins remains almost unperturbed. At the lowest value of C 
studied, that is C = 2.98, it is seen in Fig. 5d that the Sat¬ 
urn realm basin has disappeared and about 80% of the (v, y) 


4.5 An overview analysis 

It would be very informative to monitor the evolution of the 
percentages of all types of orbits as a function of the Ja¬ 
cobi constant. Fig. 6 shows such a diagram, where the verti¬ 
cal dashed black lines indicate the four critical values of C. 
We observe that at high enough values of the Jacobi constat, 
that is for C > 3.02, the collisional set is the most popu¬ 
lated family occupying more than 90% of the configuration 
space. However for C < C 2 the percentage of collisional or¬ 
bits exhibits a sharp decrease until C = C 4 , while for lower 
values of C it saturates around 12%. As the neck through Li 
opens for C < Ci the rate of escaping orbit to the Saturn 
realm increases until C = 3.06. For lower values of C the 
trend is reversed and for low enough values of the Jacobi 
constant (C < 2.98) it completely vanishes. For C < C 2 the 
second channel around L 2 also opens and it is seen that the 
percentage of orbits escaping to the exterior realm displays 
a rapid increase until C = C 4 . For lower values of C the 
slope of the curve changes and the increase is not so dra¬ 
matic. Nevertheless, when the forbidden regions disappear 
from the (v, y) plane the exterior realm basin dominates. As 
for the bounded basin it is evident that it is by far the most 
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C4 C3 C2 Ci 



Fig. 6 Evolution of the percentage of the initial conditions of each con¬ 
sidered basin as a function of the Jacobi constant. The vertical dashed 
black lines indicate the four critical values of C. 


weak type of orbits. For C > Ci the rate of bounded orbits 
has a mean value around 7.5%, in the interval C 2 < C < C 4 
its mean value drops to 4.5%, while for lower values of C 
this percentage is effectively zero because as we seen in Fig. 
5(a-d) the stability island in the left side of the Titan disap¬ 
pears. 

The color-coded planes in the configuration (x, y) space 
provide sufficient information on the phase space mixing 
however, for only a fixed value of the Jacobi constant and 
also for orbits that traverse the surface of section either di¬ 
rectly (progradely) or retrogradely. Henon back in the late 
60s (Henon, 1969), in order to overcome these limitations 
introduced a new type of plane which can provide informa¬ 
tion about areas of bounded and unbounded motion using 
the section y = i: = 0, y>0 (see also Barrio et al. (2008)). 
In other words, all the initial conditions of the orbits of the 
test particles are launched from the v-axis with v = vq, par¬ 
allel to the y-axis (y = 0). Consequently, in contrast to the 
previously discussed type of plane, only orbits with pericen- 
ters on the v-axis are included and therefore, the value of 
the Jacobi constant C can now be used as an ordinate. In 
this way, we can monitor how the variation on C influences 
the overall orbital structure the Saturn-Titan system using a 
continuous spectrum of Jacobi constant values rather than 
few discrete levels. In Fig. 7a we present the orbital struc¬ 
ture of the (v, C) plane when C e [2.995,3.042]. The four 
critical values of C are indicated by the horizontal dashed 
black lines. 


The two stability islands of prograde and retrograde mo¬ 
tion in both sides of the Titan are now visible. Here it should 
be emphasized that around the center of the Titan there is a 
horizontal lane of collisional orbits. This lane is composed 
of initial conditions of orbits started inside the radius of the 
Titan and therefore they are considered to collide with the 
smaller primary. Below the right stability island there is a 
fractal mixture of basins corresponding to escaping orbits 
to the Saturn and the exterior realm. Once more it is evident 
that for low values of the Jacobi constant the vast majority of 
the orbits escape to the exterior realm through L 2 , while for 
high enough values of C on the other hand all the initial con¬ 
ditions lead to collisional motion. Furthermore, the diagram 
shown in Fig. 7a show us how the fractality of the several 
basin boundaries strongly varies not only as a function of 
the Jacobi constant but also of the spatial variable. In partic¬ 
ular, one can observe a very interesting phenomenon. Being 
more precise we see that the fractality of the basin bound¬ 
aries, which is related to the unpredictability, migrates from 
the upper right side of the Titan for relatively high C (or in 
other words low energy) to the lower left side of the smaller 
primary for low values of the Jacobi constant (high values of 
the energy). 

In Fig. 7c we illustrate how the escape and crash times of 
orbits are distributed on the (v, C) plane. Light reddish col¬ 
ors correspond to fast escaping/collional orbits, dark blue/purple 
colors indicate large escape/collional rates, while white color 
denote stability islands of bounded motion. We note that the 
scale on the color bar is logarithmic. Inspecting the spatial 
distribution of various different ranges of the escape time, 
we are able to associate medium escape time with the stable 
manifold of a non-attracting chaotic invariant set, which is 
spread out throughout this region of the chaotic sea, while 
the largest escape time values on the other hand, are linked 
with sticky motion around the boundaries of the stability is¬ 
lands. As for the collision time we see that orbits with initial 
conditions very close to the vicinity of the center of the Ti¬ 
tan and for relatively high values of C (C > 3.03) collide 
with it almost immediately, within the first time step of the 
numerical integration. 

In order to obtain a more complete view on how the na¬ 
ture of orbits in the Saturn-Titan system is influenced by the 
energy parameter, we follow a similar numerical approach 
to that explained before but now all orbits are initiated from 
the vertical y-axis with y = yo. In particular, this time we use 
the section x = y = 0, i: > 0, launching orbits parallel to the 
v-axis. This allow us to construct again a 2D plane in which 
the y coordinate of orbits is the abscissa, while the value of 
the Jacobi constant is the ordinate. The orbital structure of 
the (y, C)-plane for the same range of values of C is shown 
in Fig. 7b. The black solid line is the limiting curve which 
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Fig. 7 Orbital structure of the (a-upper left): (x, C)-plane; and (b-upper right): (y, C)-plane; (c-lower left and d-lower right): the distribution of the 
corresponding escape and collisional times of the orbits. The color code is the same as in Fig. 4. 



Fig. 8 Evolution of the percentages of all types of orbits on the (a-left): 
The vertical dashed black lines denote the four critical values of C. 

distinguishes between regions of allowed and forbidden mo¬ 
tion and is defined as 



C)-plane and (b-right): (y, C)-plane as a function of the Jacobi constant. 


A very complicated orbital structure is reveled in the (y, C)- 
plane which however, in general terms, is very similar with 
that of the (x, C)-plane. One may observe that for C < C 4 the 
forbidden regions of motion around L 4 and L 5 completely 


fdy,C) = Q(x = 0,y) = CI2. 


( 7 ) 
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disappear. The corresponding escape and collisional times 
of the orbits are given in Fig. 7d. Looking carefully Figs. 
7(c-d) we may conclude that the smallest escape/collisonal 
periods correspond to orbits with initial conditions inside 
the escape/collison basins, while orbits initiated in the frac¬ 
tal regions of the planes or near the boundaries of stability 
islands have the highest escape rates. 

The evolution of the percentages of all types of orbits 
on the (v, C) and (y, C) planes as a function of the Jacobi 
constant are presented in Fig. 8(a-b). We observe that the 
patterns of both diagrams are very similar. At this points it 
should be very informative to briefly discuss the evolution 
of the percentages. We will analyze the patterns of Fig. 8a; 
those shown in Fig. 8b are quite similar. When C > 3.04 
all the tested initial conditions correspond to collisional or¬ 
bits. As the value of C decreases however, the rate of colli¬ 
sional orbits decreases almost linearly until Ci. In the inter¬ 
val C2 < C < Ci we observe a sudden peak at about 80%, 
while for the interval C 3 < C < C 2 the linear decrease con¬ 
tinues until C3. For all lower values of the Jacobi integral the 
percentage of collisional orbits remains constant and about 
15%. For C < Cl the neck near Li opens and the orbits may 
enter the Saturn realm. Indeed, we see that in the interval 
C3 < C < Cl the rate of escaping orbits trough Li exhibits 
fluctuations between 5% and 15%, while for C < C4 it sud¬ 
denly vanishes. The channel leading to escape to the exterior 
realm opens for C < C2, where the corresponding percent¬ 
age grows rapidly until about C4. For lower values of C it 
stabilizes around 84%. Finally the rate of the bounded basin 
increases from 0% at C = 3.038 to about 45% at Ci. Inside 
the interval C2 < C < Ci there is a sudden reduction, while 
the growth continues up to C = 3.005 (30%). For smaller 
values of C, or in other words for higher energy levels, the 
percentage of bounded orbits reduces and for C < C4 it tends 
asymptotically to zero. Therefore, one may reasonably con¬ 
clude that at low energy levels which correspond to close 
ZVCs around the Titan collisional motion dominates, while 
at high energy levels on the other hand the vast majority of 
orbits escape to the exterior region. 

5 Discussion and conclusions 

The planar version of the circular restricted three-body prob¬ 
lem was used in this contribution in order to numerically in¬ 
vestigate the orbital dynamics of a small body (spacecraft, 
comet or asteroid) under the influence of the potential of the 
Saturn-Titan system. All the initial conditions of orbits were 
initiated in the vicinity of the Titan which was our scattering 
region. We managed to determine four basins corresponding 
to: (i) bounded orbits around the Titan, (ii) escaping orbits 
to the Saturn realm through Li, (iii) escaping orbit to the 
exterior realm through L 2 , (iv) leaking orbits due to colli¬ 
sions with the surface of the Titan. The orbital structure of 


the basins was monitored by varying the value of the Jacobi 
constant. As far as we know, this is the first time that the or¬ 
bital content in the vicinity of the Titan is explored in such 
a detailed and systematic way. 

We believe that our results could be useful in space mis¬ 
sion design. Nevertheless, we have to stress out that in real 
space missions there are several types of plane perturbations 
mainly due to the presence of other heavy celestial bodies 
(e.g., planets). Near the vicinity of Titan (which was the 
scattering region in our work) however, all these perturba¬ 
tions are extremely weak and therefore negligible, so we 
may claim that the presented numerical results are struc¬ 
turally stable against out of plane perturbations. 

For several values of the Jacobi integral we defined dense, 
uniform grids of 1024 x 1024 initial conditions (vo,yo) reg¬ 
ularly distributed in the area allowed by the Jacobi constant 
on the configuration (v, y) space and then we forwarded in¬ 
tegrated them. For the numerical integration of the orbits in 
every grid, we needed about between 1 minute and 5 hours 
of CPU time on a Pentium Dual-Core 2.2 GHz PC, depend¬ 
ing on the escape, collisional and bounded rates of orbits in 
each case. For each initial condition, the maximum time of 
the numerical integration was set to be equal to 5000 dtu 
however, when a particle escaped or collided the numerical 
integration was effectively ended and proceeded to the next 
available initial condition. 

Our work contains quantitative information regarding the 
orbital dynamics in the Saturn-Titan. The main numerical 
results can be summarized as follows: 

1. In all examined cases, regions of initial conditions lead¬ 
ing to escape in a given direction or to collision with 
the smaller primary were found to exist. The several es¬ 
cape/collision basins are very intricately interwoven and 
they appear either as thin elongated spiral bands or as 
well-defined broad regions. 

2. It was found that for high values of the Jacobi constant, 
or in other words for low energy levels, the configuration 
space is dominated by collisional orbits, while some sta¬ 
bility islands are also present. 

3. We observed that for low values of the Jacobi constant 
which correspond to high values of the total orbital en¬ 
ergy the vast majority of the orbits escape to the exterior 
realm, while the rest of the (x, y) plane is occupied by 
collisional orbits. 

4. It was detected that the average collisional time increases 
with increasing C, while on the other hand, the average 
escape time, in general terms, diminishes with decreas¬ 
ing C. 

5. Our calculations revealed that inside the basins the short¬ 
est escape and collision rates of the orbits had been mea¬ 
sured. On the contrary, the longest escape/collision rates 
correspond to initial conditions located near the bound- 
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aries between the basins or in the vicinity of the stability 
islands. 

Taking into account the detailed and systematic analy¬ 
sis we may conclude that out numerical task has been suc¬ 
cessfully completed. We hope that the present results to be 
useful in the active field of orbital dynamics in the restricted 
three-body problem. Since our outcomes are encouraging, it 
is in our future plans to expand our investigation into three 
dimensions and also in other planetary systems. 
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